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In his Comment |arXiv: 1501.02126 (2015)] on our recent paper [Phys. Rev. Lett., v. 113, 254101 
(2014)], Pikovsky compares two methods for defining the “phase” of a stochastic oscillator. We 
reply to his Comment by showing that neither method can unambiguously identify a unique system 
of isochrons, when multiple oscillations coexist in the same system. 


In his comment [T] on our paper [5], Pikovsky con¬ 
trasts two definitions for the phase of a stochastic oscil¬ 
lator by way of an analytically solvable model system. In 
[3] the phase is defined in terms of a system of isochrons 
Eg, analogous to Poincare sections, with the property 
that for any initial condition on one isochron E^, the 
mean first passage time (MFPT) to a second isochron 
Ef,, & > a, will depend only on the phase difference b — a. 
In our approach [3] the phase is defined as the complex 
argument of the slowest decaying eigenfunction of the 
backward Kolmogorov operator, provided the first non¬ 
trivial eigenvalue is complex and is well separated from 
the next slowest decaying eigenvalue. Pikovsky argues 
that our eigenfunction approach does not properly work 
in all situations and proposes the following example to 
demonstrate this. Consider two independent phase-like 
variables, each taking values in [0, 27r), that obey 

01 = Wi -I-CTlCl, &2=^2+<^2^2, (1) 


This quantity satisfies an equation involving the same 
backwards operator as that identifying the isochrons in 
[2], namely 


dr 2 


( 2 ) 


with boundary condition T(a, 5) — 0 as a —>■ 6 [4] . Clearly 
the solution is r(a, b) = {b— a)/uj. Therefore the surfaces 
61 = const, provide a system of MFPT isochrons for the 
system Q, as described in [3]. However, so do the sur¬ 
faces 02 = const. Moreover, for any nontrivial pair of 
integers (n,m), the surfaces ifn.m = const, form another 
system of MFPT isochrons, where we define 


if 


n,m 


nOi +17102 
71 + m 


(3) 


This is easily seen, since ifn,m obeys a stochastic differ¬ 
ential equation (SDE) of the same form as 0i and 02, 


where = 2Sij6{t — t'). The eigenvalues of 

the adjoint Fokker-Planck operator are Xn,m = *(nwi -I- 
mu) 2 ) — (ntJi)^ — {ma 2 )^. The “slowest decaying mode” 
corresponds either to Ai^o (for cti < CT 2 ) or to Ao,i (for 
a I > cr 2 )- In the first case, 0i could be interpreted as the 
primary phase variable; in the second case, 02 could be. 
(We adopt the notation 0i,02, rather than 0, (f, so as not 
to prejudge the roles of the variables.) If tJi « tT 2 , the 
system is not “robustly oscillatory” in the sense of [2], 
and neither variable is clearly identified as the primary 
“phase”. 

Thus, as Pikovsky’s example illustrates, the spectral 
method does not unambiguously identify a unique phase 
variable, when multiple oscillations coexist with similar 
coherence times. Pikovsky asserts in his comment that 
“the approach of Ref. 0 yields here the proper phase 0 — 
const.” However, as we demonstrate below, the MFPT 
method introduced in [5] necessarily exhibits the same 
ambiguity under the same circumstances, at least for this 
simple example. 

Let ijj S [0, 27r) satisfy ip = uj + a^, where ^ is white 
Gaussian noise and a; > 0. Let r(a, b) be the MFPT of 
the system starting at '0 = a to arrive a,t if = b > a. 


■ nuji + mu!2 Wn^aT+rri^a^ 

Vn,m = -^- h - , -6- 


n + m 


(4) 


Therefore, there is a countably infinite collection of sur¬ 
faces satisfying the MFPT property for the system he 
describes. It is difficult to see how the MFPT approach 
would identify a unique system of isochrons, without be¬ 
ing supplemented by additional criteria. 

We have omitted the radial variable from our eq. Q. 
In equations (1-3) of Pikovsky’s comment we note that 
the radial variable is entirely uncoupled from the two 
phase variables. The physical motivation for the example 
is a noisy limit cycle tracing an orbit in three dimensional 
space, rotating simultaneously in both angles describing 
the points on a torus. However, neither the SDE nor the 
Eokker-Planck equation analyzed in the note correspond 
to this physical system. What is missing is the interac¬ 
tion of one of the phase variables with the radial variable 
in the SDE. By decoupling them in the equations given, 
Pikovsky has made the system symmetric with respect 
to exchange of 0i O 02 (equivalently, 0 ^ (f). The sys¬ 
tem lacking this symmetry is more difficult to analyze. 
Certainly one could construct a 3D system in which the 
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FIG. 1: Multirhythmic behavior in Izhikevich’s low-threshold persistent sodium plus potassium model, with parameters giving 
subcritical Hopf bifurcation near I ~ 43/rA/cm^. For Ntot = 100 discrete, randomly gated potassium channels, and 7 = 
40pA/cm^, trajectories show a mixture of spiking and subthreshold oscillations, at two distinct frequencies. (A) Voltage 
(mV) and (B) N , number of open potassium channels are plotted against time (msec). (C) The trajectory in the {V, N) plane 
shows coexistence of a large and a small oscillation. For this value of 7, the eigenvalue spectrum computed in [2] shows two 
slowly decaying complex modes with similar negative real parts. Trajectory generated via an exact simulation algorithm [^. 


asymptotic phase obtained from the adjoint eigenfunc¬ 
tions appears ambiguous; however in any such system we 
suspect that the construction based on MFPT isochrons 
will suffer from the same ambiguity. 

Indeed, such ambiguity arises naturally in the case 
of multirhythmic or mixed-mode oscillations. In [5] we 
analyzed the eigenvalue spectrum of Izhikevich’s “low- 
threshold persistent sodium plus potassium” model, with 
parameters giving a subcritical Andronov-Hopf bifurca¬ 
tion for injected current 7 « 43/rA/cm^ ([3, Figs. 6.16 
and 4.1b). Well above the bifurcation point (at I = 
60^A/cm^) the eigenvalues follow a nearly parabolic 
spectrum, as expected for a robustly oscillatory system 
with weak noise. Just below the bifurcation point (at 
7 = 40/rA/cm^), channel noise induces switching between 
spiking and subthreshold oscillations (Fig.[^. The eigen¬ 


value spectrum shows two distinct slowly decaying modes 
with similar decay rates, reflecting the coexistence of dis¬ 
tinct oscillatory processes, each with its own typical fre¬ 
quency. It is not clear how the MFPT based isochron 
construction would perform in this setting; in our opin¬ 
ion both approaches merit further development. 
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